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Abstract 

We have studied the molecular scale hydrophobicity of an apolar solute, argon, confined between 
hydrophobic planar surfaces with different confinement widths. Specihcally, we find that the hy¬ 
drophobicity exhibits a non-monotonic behavior with confinement width. While hydrophobicity 
is usually large compared to bulk value, we find a narrow range of confinement width where the 
hydrophobicity displays similar values as in bulk water. Furthermore, we develop a simple model 
taking into account the entropic changes in nanoconfined geometry, which enables us to calculate 
potential of mean force between solutes as the conditions change from bulk to different degrees 
of planar nanoconfinement. Our results are important in understanding nanoconfinement induced 
stability of apolar polymers, solubility of gases, and may help design better systems for Enhanced 
Oil Recovery. 
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I. INTRODUCTION 


Confined systems are prevalent in natnre. Examples of confinements inclnde intracellnlar 
water and biomolecules to gases and liquids trapped in porous rock sediments One can 

expect that when conhnement length scales approach the inherent characteristic molecular 
length scales, such as various correlation lengths, the properties of a system may deviate 
from bulk [1]. 

It is widely believed that besides the interaction of solvent with solutes, thermodynamic 
states and dynamics of solvent, water in case of aqueous systems, are important in dictating 
the behavior of solubility of solutes. Nanoconhnement induced changes in the properties of 
water is quite extensively studied in recent times [5HTB]. It has been shown that when water 
is conhned to length scales of few molecular diameters, water exhibit distinct structural and 
dynamical changes as compared to bulk water im. Koga et. ah have studied liquid-solid 
phase transitions in water conhned between hydrophobic surfaces and have shown that for 
appropriate dimensions of the hat hydrophobic conhnements, water can undergo a hrst- 
order phase transition from a bilayer liquid to a bilayer hexagonal ice UU- Further studies 
of water in nanoconhnements have suggested a range of exotic structural phases such as 
n-fold crystal (n=2,4,6, and 12) in bilayer water [TOj. Bai et. ah have found a variety of 
polymorphic and polyamorphic transitions in bilayer water conhned between hydrophobic 
surfaces [18]. Recently, Han et. ah have studied liquid-solid phase transitions in bilayer 
water conhned between planar hydrophobic surfaces m- They hnd that nature of liquid- 
solid phase transition changes from a hrst-order transition at low densities to a continuous 
liquid-solid transition at high densities. Besides the structural changes in nanoconhnement, 
it has been shown that the physical property such as temperature of maximum density 
(TMD) of water shifts to lower temperature and lower density in the temperature-density 
plane as compared to bulk water [1]. Gallo et. ah have studied the ehect of nanoconhne¬ 
ment on the dynamics and thermodynamics of water and found that the water conhned in 
nanopores exhibits diherent dynamical regimes depending on the distance of water from the 
surface [9|. The ehects of morphology and charges of conhning surfaces have been explored 
in Ref. [HI Il9| . 

Aqueous solubility and hydrophobicity of solutes is an active area of research due to 
its implication in a large variety of physical phenomena including protein stability and 
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folding [201 - I26] . hydrate formation [27], and enhanced oil recovery [3|. Studies of solubility 
of hydrocarbon solutes in water suggests that they exhibit anomalous solubility in water pSh 
150] . For example methane solubility in water shows a minimum around 320 K and increases 
on both sides of this temperature. Low temperature increase of solubility of apolar solutes 
is implicated in cold unfolding of proteins and apolar polymers [25] . In more recent works, 
Mallamace et. al. have explored the influence of water on protein properties including 
protein folding and dynamic transition in proteins at low temperatures na Eg. Recent 
computational studies have investigated the effect of nanoconhnement on the phase behavior 
of oil-water mixture in an effort to improve the existing enhanced oil recovery technologies [3] . 
Studies of solubility and ordered phase formation of gases in solid ice phases have attracted 
wide attention due to its importance in hydrocarbon processing and sustainable energy 
production [271 ED ED El] ■ 

Although a large body of literature exists on the solubility and hydrophobicity of apolar 
solutes in bulk water, the effect of nanoconhnement on the solubility and hydrophobicity 
remains poorly understood. Since both thermodynamic and dynamic properties of water 
may vary in nanoconhnement, one may expect that the hydrophobicity may also be diherent 
from the bulk . Water in biological systems is under the conditions which diher greatly from 
the bulk hence the bulk behavior of water may not be relevant under those conditions. 
Therefore, a good understanding of the behavior of hydrophobicity of any molecule would 
only be possible if a detailed study of the these conditions are explored. 

In this article, we study molecular scale hydrophobicity of an apolar solute, argon, in pla¬ 
nar hydrophobic nanoconhnement of varying length scales. In Method section, we describe 
the simulation methods, in Results section, we discuss the results and hnally we conclude 
with Summary and Discussion section. 


II. SYSTEM AND METHOD 

We performed molecular dynamics simulations to calculate potential of mean force 
(PMF) [HE] between argon atoms dissolved in TIP3P (transferable interaction potential 
three points) [35l 136] water-like solvent conhned between two structured planar hydrophobic 
surfaces separated by a distance L^. The atoms on the surfaces were arranged in a hexagonal 
closed packing with a distance as = 0.339 nm and they do not interact with each other. 
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The positions of the surface atoms were restrained to their respective mean positions by a 
harmonic potential with a force constant 10^ kJ/mol. In order to mimic hydrophobic sur¬ 
faces, carbon atoms were chosen to be surface atoms. The interaction between the surface 
atoms and the oxygen of the water molecule is modeled using 6 — 12 Lennard-Jones (LJ) 
potential U (r) 

U{r)=Ae[{^r-{^f] ( 1 ) 

^SO rso 

where rso is the distance between the surface atom and oxygen of water molecule, and 
f = (0‘s + o'ow)/2 and e = yTshow- o'ow and eow are the parameters of the Lennard-Jones 
interaction between oxygen atoms of water molecules [351 136] . 

We performed simulations of conhned water and solutes between hydrophobic surfaces in 
NVT-ensemble with effective density of the system, p = 1.00 g/cm^. The effective density 
of the system is obtained by calculating the effective confinement width available for water 
molecules for a given center-to-center distance, L^, between the surface atoms. The effective 
conhnement width, available to water molecules is calculated as 




(O‘ow + O'5o) 
^ 2 


( 2 ) 


The equations of motion were integrated with a time step of 0.001 ps and velocity rescaling 
was used to attain constant temperature and Berendsen barostat for constant pressure in 
Gromacs 4.5 EH. Periodic boundary conditions were applied in XY-directions. We use 
constrained molecular dynamics method for the calculation of PMF for seven different ef¬ 
fective conhnement widths ^ = 0.6, 0.8,1.0,1.2,1.4,1.6, and 1.8 nm. Images of some of the 
representative conhgurations studied here are shown in Fig. [^a). 

After the system was equilibrated for 1 ns at T = 300 K for respective conhnements, 
we constrain the argon atoms at hxed values of distances with a harmonic potential with a 
spring constant kc = 1000 kJ/mol. After the equilibration step, we ran the simulations for 
additional 1 ns for each constraining distances d between 0.26 and 1.0 nm at an interval of 
0.02 nm. Potential of mean force was corrected for volume entropy at T= 300 K. For the 
range of conhnement widths studied here, we hnd strong layering of water near the surfaces 
as suggested by the transverse density prohle of water along the conhnement direction z (see 
Figure [^b)). 

While the ehective density of the system remains the same, the lateral pressure along 
the periodic directions, Pxy, varies for diherent conhnement widths as shown in Fig. [^ P^y 
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monotonically increases with effective confinement width 


III. RESULTS 

POTENTIAL OF MEAN FORCE AND SECOND VIRIAL COEFFICIENT 

To quantify hydrophobicity, we first calculate PMF, ta(r), for different conhnement widths 
using umbrella sampling [SH]- To use umbrella sampling, the distance between the argon 
atoms is chosen as the reaction co-ordinate. Within this scheme, a constraining potential is 
added to the Hamiltonian. The modified Hamiltonian of the system with the constraining 
potential with a perturbation parameter A can be written as 

+ (3) 

where kc is the spring constant of the constraining potential and tq is the constraining 
distance. The free energy change AG between two arbitrary points along the reaction 
coordinate is given by 

Hence the free energy difference between two states along the reaction coordinate is integral 
of the mean force between the state points. The PMF, w{r), for a given state point along 
the reaction coordinate with respect to a reference state is 

w{r) = — f < F{rQ) > drQ + 2kBT\n{r) + C. (5) 

where, < Flr'^) > is the mean force for a given constraining distance Tq, Uref is the position 
along the reaction coordinate for the reference state, and the second term in the above 
expression is the volume entropy correction. We choose the reference state to be = 1 nm 
where w{r.^ei) = 0. 

In Fig.|^ a) and (b), we show PMF, w{r), as a function of distance r between argon atoms 
for all the conhnement widths. For a comparison, we also show the PMF for argons in bulk 
water at P = 1 atm and T = 300 K. w{r) for bulk system exhibits two prominent minima, 
one at r 0.35 nm and another at r ^ 0.68 nm respectively. While the position of the hrst 
minimum of w{r) remains unchanged in conhnement, the position of the second minimum 
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decreased to r 0.64 nm for the smallest Moreover, while the hrst minimum becomes 
deeper monotonically with decreasing conhnement width, the second minimum shows a 
non-monotonic behavior with ^ suggesting a non-trivial solvation structure dependence with 
conhnement width 

In Fig-I^c ) and (d), we show the radial distribution function between argons as calculated 
from gAr-Ar{r) = for all the conhnement widths studied here. It is clear that solvation 

structure of argons in water in nanoconhnement exhibits a non-trivial dependence on the 
conhnement width of conhning surfaces. Since, this non-trivial behavior of PMF or gAr-Ar{f') 
for the second shell makes it harder to interpret the hydrophobicity directly by looking at 
them, therefore to quantify hydrophobicity, we next calculate the second Virial coefficient. 

The second Virial coefficient B 2 for a mono-component system is given by 

1 r°° 

B2 = -a - l)dF (6) 

4 Jo 

where (3 = and dT is the volume element corresponding to the separation r. w{r) is 
related to the solute-solute pair correlation function h{r) as 


w{r) =—kBTln[h{r) + 1] (7) 

and hence 

B 2 = -i/ Kr)dT ( 8 ) 

A positive value of B 2 indicates an ehective repulsive interaction between the solutes and 
hence larger solubility and a negative value suggest an ehective attractive interaction and 
hence smaller solubility. The volume element, dF = A{r,^)dr, in the planar conhnement 
depends on both r and ^ and can be derived by hrst deriving the average surface area A{r, 
corresponding to solute separation r. Let’s assume that the surfaces are represented by two 
planes; one at z = 0 and another at z = ^. Assuming that the hrst particle can assume 
any value of between 0 and we can next write the average area, A(r;^), traversed by 
the radial vector joining two particles as (a detailed derivation of A(r;^) is given in the 
Supplementary Information): 


A%r;0 


( 


< 


V 


47rr^[l 

27rre(l - (|)2) +47rr2(l 


- -1 

if r < ^/2 

- —) 
2rl 

if ^/2 < r < ^ 

27rr^ 

if r > ^ 


(9) 
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The above expression suggests that the accessible area to arrange a pair of particle of fixed 
separation r decreases with r as compared to bulk. 

In Fig. 1^ we show the second Virial coefficient B 2 as a function of The values of 
B 2 are signihcantly smaller compared to the bulk value of B 2 {—17A^) at P = 1 atm and 
T = 300 K except for .^ = 1.2 nm, suggesting that the hydrophobicity of argon increases 
in nanoconhnement. Moreover, B 2 exhibits a nonmonotonic dependence on B 2 hrst 
increases with increasing ^ for ^ = 0.6 — 1.2 nm and then it falls off sharply with a weak 
dependence on ^ for .^ = 1.4 — 1.80 nm. 

Since the simulations are done in NVT-ensemble, the pressure for different values of ^ are 
different even when the effective density remains the same. One may argue that the change 
in hydrophobicity is just due to pressure increase as the conhnement width decreases. To 
test this we next calculated vlaues of B 2 for argon for two different pressures P = 1 and 
500 atm in bulk water. The pressures were chosen such that they cover the pressure ranges 
for the nanoconhned system. In Fig. (a), we show the PMF for P=1 and 500 atm for the 
bulk system at T = 300 K. The second Virial coefficient B 2 for P = 1, 500 atm are —17 and 
53 A° respectively. We hnd that the bulk values of B 2 differ greatly from the nanoconhned 
values except for ^ = 1.2 nm, for which the system shows a comparable value of B 2 . To 
this end, it clear that the magnitude of hydrophobicity in nanoconhnement is drastically 
diherent from the bulk. Moreover, our results indicate subtle changes in solvation and 
ehective hydrophobicity upon nanoconhnement. Indeed, small changes in the length scale 
of nanoconhning regions may lead to large changes in hydrophobicity. 


FROM BULK PMF TO CONFINED PMF 


Nanoconhnement results in changes in arranging a solute particle around a reference 
particle (see Supplementary Information) -namely a decrease in the volumetric arrangement, 
we can dehne the change in entropy AS'i(r;.^) in nanoconhnement over a bulk system as: 


A5i(r;0 


A:Bln(l-^) if0<r<e/2 
< fcBln[(l-|:) + |:(l-(|)2)] if e/2 < r < ^ 

fcslnd:) ife<’^<’^o 


( 10 ) 
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where ro is the separation between solutes such that hnir_>r.o'?i'(?’) = 0. Note that ro is 
the distance beyond which the system loses two-point correlations. In our case, we assume 
ro = 1 nm. The validity of above expression ranges for the values of r over which the solute 
particles are correlated beyond which the two point entropy would be zero for both bulk and 
conhned system. Since the entropy decreases as r is increased, we argue that the system with 
solute nano confinement will try to increase the configurations pertaining to larger entropy in 
order to minimize free energy. Hence, we can assume an additional thermodynamic driving 
potential —TAS'i(r;^) acting on a pair of particles. This entropic penalty will make small 
distances between solutes more favorable as they correspond to higher entropy. Taking this 
into account, we can write the modihed PMF in nanoconhnement w*{r]^) as 

w*(r;0 = Wbuik(r)-TA^i(r) (11) 

The w*{r;^) is then normalized such that w*{r = Inm;^) = 0. In Fig. we compare PMF 
computated from simulations and PMF from Eq. [T^for confinement widths ^ = 0.60 nm, 
^ = 0.80 nm, ^ = 1.60 nm, and, ^ = 1.80 nm. We hnd that the PMF predicted by theory 
does reasonably well for smaller ^ but it deviates from the PMF calculated from simulations 
for larger The deviation of theory could presumably results from absence of enthalpic 
consideration. We show the values of B 2 calculated from w*{r]^) for different values of 
conhnement widths in Fig. along with the computed values of B 2 from simulations. 

We hnd that simple argument of entropy decrease of conhgurations with different sep¬ 
aration of a pair of solute atoms in conhnement works well with planar geometry studied 
here. Indeed, a similar argument can be made for diherent geometries of conhnement such 
as cylindrical nanoconhnement. 


SOLVATION STRUCTURE 

Next, we investigate the solvation structure in order to hnd a clue to the anomalous 
hydrophobicity behavior for ^ = 1.2 nm by investigating the radial distribution of water 
and argon. For water, we calcualte the oxygen-oxygen lateral radial distribution function, 
gooi^xy) as a function their lateral distance r^y In Fig. (a), we show gooir^y) for all 
the conhnement widths We hnd that, water tends to order laterally with the decrease 
of conhnement width as suggested by increased hrst and second peaks in 5'oo(^xy) upon 


decrease of In Fig. [^b), we show the radial distribution function between oxygen of 
water and argon, go-Ari^) calculated from long simulations of water-argon system (4 ns for 
each conhnement width). While the structure of the hrst solvation shell (r 0.36 nm) does 
not show any appreciable change with ^ except the decreasing value of the hrst peak, the 
second shell becomes much wider and the corresponding peak moves to slightly larger values 
of r. Combining the results on PMF between argon and solvation structure, it seems that 
non-monotonic dependence of second shell on ^ may explain the increased value of B 2 for 
^ = 1.2 nm and will require further investigation. 

IV. SUMMARY AND DISCUSSION 

We have studied molecular scale hydrophobicity of a small apolar solute in nanoconhne- 
ment by calculating the second Virial coefficients from the potential of mean force between 
apolar solutes in water conhned between hydrophobic surfaces with different conhnement 
widths. We hnd that: (i) hydrophobicity of apolar solutes usually increases as the con¬ 
hnement width decreases at constant transverse pressure, (ii) hydrophobicity exhibits an 
anomalous region, where the hydrophobicity becomes similar to the values in bulk wa¬ 
ter. Furthermore, we hnd that this anomalous region of conhnement widths correspond 
to changes in second and third solvation shell. We develop a simple entropic theory to 
hnd ehective potential of mean force in nanoconhnement. The predicted PMF from theory 
works reasonably well for smaller conhnement widths and allows to predict the ehective 
hydrophobicity in nanoconhnement when the potential of mean force in bulk is known. 

Hydrophobicity plays important roles in many physical systems including the folding of 
polypeptides. We hypothesize that large changes in hydrophobicity with very small changes 
in length scale of nanoconhning regions could be very important in maintaining the rate of 
folding in biological systems such as chaperone proteins in a chemically non-specihc way. 
A very good example of this is GroEL-GroES complex in bacteria, which is induced upon 
temperature shocks [SH]- To counteract the increase of temperature and hence unfolding of 
polypeptides, nature has evolved a complex machinery. An unfolded polypeptide is directed 
into a chaperon complex such as GroEL-GroES and stays inside the cage formed by the 
protein complex where it folds much faster than it would in the bulk region [58]. We argue 
based on our results that when and unfolded or partially folded polypeptide enters this 
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complex, a large part of the folding may occur in the neck-region before it reaches the cage. 
Moreover, we hypothesize that the modulating diameter of this chaperon complex would 
naturally be helpful in modulating the rate of folding in a chemically non-specihc way. 

In summary, our results are important in understanding nanoconhnement induced stabil¬ 
ity of apolar polymers, solubility of gases and may help design better systems for enhanced 
oil recovery. While our work has explored the planar hydrophobic conhnement in details, the 
geometry and charge distribution of the conhning surfaces may affect the hydrophobicity. 
Furthermore, effect of solute size consideration is also a determining factor for hydrophobic¬ 
ity which we will explore in future work. 
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FIG. 1. (color online) (a) Schematic of nonoconfined system, (b) Potential of mean force (PMF) 
between two argon atoms as a function of distance r between argon atoms for different confinement 
width 
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FIG. 2. Lateral pressure P^y as a function of effective confinement width ^ for the fixed density 
p = l.Og/cm^. Pxy decreases monotonically with confinement width. 
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FIG. 3. (color online) Potential of mean force between argon in nanoconfiment for different confine¬ 
ment widths (a) ^ = 0.6, 0.8,1.0, and (b) ^ = 1.2,1.4,1.6,1.8 nm respectively. For a comparision, we 
also show PMF for argon in bulk water at T = 300 K and P = 1 atm. Argon-argon radial distribu¬ 
tion function (7Ar-Ar(?') calcuated using PMF w{r) for confinement widths (c) = 0.6,0.8,1.0 nm, 

and (d) ^ = 1.2,1.4,1.6,1.8 nm. 
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FIG. 4. (color online) Second Virial coefficient B 2 between argons as a function confinement width 
Also shown is the curve obtained using entropic consideration in confinement (see Section V). 
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FIG. 5. (color online) (a) Potential of mean force between argon, w{r), for bulk system for two 
different pressures P = 1, and 500 atm and temperature T = 300 K. (b) Second Virial coefficient 
B 2 for two pressures P = 1, and 500 atm and temperature T = 300 K. 
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FIG. 7. (color online) Comparision of potential of mean force calculated from simulation with 
theory for confinement widths (a) ^ = 0.60 nm, (b) ^ = 0.80 nm, (c) ^ = 1.60 nm, and (d) 
^ = 1.80 nm. While the PMF predicted from theory does reasonably well for smaller it deviates 
from the PMF calculated from simulations for larger 


19 








































1.8 

1.6 

1.4 

1-2 

o 

“ 1.0 
0.8 
0.6 


0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

r fnml 

xy 



(a) 



n ^=0.60nm 
3 ^=0.80nm 
^=1.00nm 
‘^=1.20nm 
<^=1.40nm 
''^= 1 . 60 nm 
^=1.80nm 


FIG. 8. (a) Lateral oxygen-oxygen radial distribution function, 5oo(?’xy) of water for various con¬ 
finement width, (b) Oxygen-Argon radial distribution function, go-x^(r), for all the confinement 
widths studied here. 
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Supplementary Information: 


Effective surface area traced by a particle at a distance r from another particle 

Let’s assume that two infinite planar surfaces are located at z = 0 and z = ^ respectively. 
We are interested in the average surface area traced by a particle at a fixed distance r from 
the reference particle. It is easy to see that the average surface area for a fixed r depends 
on both r and When the reference particle is closer to the surface the area is smaller 
compared to when the reference particle is sitting close to the center of the surfaces. Three 



Z 






Case (iii) 


FIG. SI. Schematics of the regions of r discussed above for a confinement width 

different regimes of r - (i) 0 < r < ^/2, (ii) ^/2 < r < ^, and (iii) r > ^ can be noted. 
Furthermore, from the symmetry of the system, it is sufficient to calculate the average area 
traced by the second particle for different positions of the reference particle only along the z 
direction. Let’s assume that the reference particle is sitting at a distance from the surface 
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at z = 0. In the following, we will derive the expression for the average surface area traced 
by a second particle placed at a hxed distance r from the reference particle. 

Case (i) ; 0 < r < ^/2 

The average area is given by 

1 

Mr-,0 =c Mr,z-,Odz 
s Jo 

where A{r, z-,0 is the surface area traced by the second particle when the reference particle 
is located at a distance z along the z-direction. For this case, we can split A(r, z;0 ii^fo f^o 
terms, (a) when z>r or ^ — z>r and (h) z < r or ^ — z < r (see Fig. [Slj ). Hence, the 
total area A{r,z]0 can be written as 

nr 

H(r, z; = 47rr^(^ — 2r) + 2 / 27r(r + z)dz = 47rr^^(l-) 

Jo 2,^ 

And hence, the average area A(r, .^) can be written as 

Air-,0 = 47rr2(l - 


Case (ii) ; i/2 < r < ^ 

In this case, we can show that the total area A(r, .^) can be written as a combination of 


two terms when [a) z < i — r, and (b) ^ — r < 2 : < i/2 (see Fig. SI), leading to 


^ C 2 

27rr(r + z)dz + 27rr^(r - i/2)] = 27rr^(r - |) + 7rr^^(l - 
And hence the average area, A(r;^) is given by 


A(r;0 = 47rr2(l “ ^) + 27rr^(l - 


Case (iii) : r > i 


In this case the total area Air, z-,0 is independent of the location 2 ; of reference particle 


and is 27rri and hence the average area A(r;.^) is given by (see Fig. SI) 


A(r;0 = ‘2nri 
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In Figure S2 (a), we show the effective average area as a function of r for all the 

confinement widths studied here. Also, in Figure [S^b), we show entropy change AS'i(r, .^) 
as a function of r. Entropy change, AS'i(r;^), over the bulk system is defined as: 


ASi(r;0 = 


( 12 ) 


where ks is the Boltzmann constant. AS'i(r;^) represents the change in entropy of vol¬ 
umetric arrangement of a two-particle system in confinement over the bulk situation. As 
expected, the entropy change for small r in confinement is smaller compared to large r which 
gets increasingly larger as the confinement width is reduced. 
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FIG. S2. (a) Effective average area traced by a particle for a fixed distance r from a reference 
particle, and (b) Entropy change AS'i(r;^) as a function of r for different confinement widths. 
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